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Abstract 

Frangais. 

La methode d'imagerie passive en sismologie a ete developpee recemment 
en vue d'imager la croute terrestre a partir d'enregistrements du bruit sis- 
mique. Elle repose sur le calcul des fonctions de correlation de ce bruit. 
Nous donnons dans cet article des formules explicites pour cette correlation 
dans le regime "semi-classique" . Pour cela, nous definissons le spectre 
de puissance d'un champ aleatoire comme I'esperance de sa mesure de 
Wigner, ce qui permet d'utiliser un calcul dans I'espace des phases : le 
calcul pseudo-differentiel et la theorie des "rays". Nous obtenons ainsi 
une formule pour la correlation du bruit sismique dans le regime "semi- 
classique" avec une source de bruit qui pent etre localisee et non homogene. 

Nous montrons ensuite comment I'utilisation des ondes guidees de sur- 
face permet d'imager la croute terrestre. 

Mots des : Imagerie passive ; semi-classique ; ondes de surfaces. 

English. 

The method of passive imaging in seismology has been developped re- 
cently in order to image the earth crust from recordings of the seismic 
noise. This method is founded on the computation of correlations of the 
seismic noise. In this paper, we give an explicit formula for this correlation 
in the "semi-classical" regime. In order to do that, we define the power 
spectrum of a random field as the ensemble average of its Wigner measure, 
this allows phase-space computations: the pseudo-differential calculus and 
the ray theory. This way, we get a formula for the correlation of the seismic 
noise in the semi-classcial regime with a source noise which can be localized 
and non homogeneous. After that, we show how the use of surface guided 
waves allows to image the earth crust. 
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Introduction 



Correlations of the noisy wave fields is used as a new tool in seismic imaging and 
monitoring, starting from the pioneering work of Campillo and Paul [3] (similar 
tools have been used in helio-seismology [11]) and followed by many works [TJ El 
[m [HI HH [201 [211 1271 [28] . See also the review paper [H] . It has also been used in 
the monitoring of the deformations of volcanoes [2]. Because it is a very powerful 
method and, hopefully, in order to make it more efficient, it is quite challenging to 
give mathematical supports to this method, now called "passive imaging". This 
has been done in a rather great generality in [6l [5] using semi-classical analysis 
(see also [ISl [13 [B [H]). 

Exact formulas for the correlations of the fields are known if the source noise 
is homogeneous (a white noise). This assumption is not satisfied in applications. 
It is therefore desirable to get formulate valid for more general source noises, 
in particular if the source noise is localized in some part of the domain. This 
turns out to be possible in the so-called semi-classical regime where the wave- 
lengths are negligible with respect to the size of the propagation domain. The 
field correlation admits a general expression in terms of the Green's function and 
the source correlation (Equation (^). The idea is to find the asymptotics of this 
expression in the semi-classical regime. 

I will present in this paper approximate formulas which are valid in the range 
of high frequency wave propagation and for which the source noise is localized 
in some part of the domain of propagation. The correlation is explicitly given 
in term of the decomposition of the Green's function as a sum over rays and 
the (phase-space) power spectrum of the source noise. I can use ray theory if I 
assume that the source noise has a short correlation distance of the same order 
of magnitude than the wavelengths. This calculus can be presented in a very 
geometric way using rays propagation as well as a re-interpretation of the source 
correlation in terms of the phase space power spectra. I use the calculus of 
pseudo-differential operators in a very essential way. I will not reproduce the 
mathematical arguments which are presented in my paper [5], but I will try not 
only to give explicit formulas, but also to present the main ideas and tools. 

Here is a more precise description of the content: the goal is to get the formula 
given in Theorem 14.11 which gives the modification of the correlation of the seismic 
noise induced by the non-homogeneity of the source noise. The modification is 
given in terms of the power spectrum of the source noise, the attenuation and 
the ray dynamics associated to the deterministic wave equation. 

I first give a review of the pseudo-differential calculus (section [1]): this allows 
to put the basic terminology of rays dynamics and to define power spectra of 
arbitrary random fields (section [2]). 

I then introduce the simplest mathematical model where the source noise is 
simply the right-handside of the wave equation (section [3]) and I present our main 
formula in section [H The interest of the result depends of the relation between 
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2 time scales discussed in section |5l the Ehrenfest time given in terms of the 
Lyapounov exponent and the attenuation time. 

How to use all of this in imaging problems? I do that (section |6]) in the case 
of seismology using the effective wave equation for the guided surface waves. The 
final problem turns out to be an inverse spectral problem whose mathematical 
solution is known. 

Finally, I discuss in section [7] a related issue, namely the calculus of the 
correlations of plane waves scattered by an obstacle or an inhomogeneity viewed 
as random waves: the direction of the waves is supposed to be random and 
uniform. This way, I show that the result of [I9] is completely general. 

1 A short review of the pseudo-differential cal- 
culus and Wigner measures 

For the mathematics of pseudo-differential operators, see P, [T21 [131 121] • 

The pseudo-differential operators (\I^DO's) were introduced in the sixties by 
Calderon, Zygmund, Nirenberg, Hormander and others as a tool in the study of 
linear partial differential equations with non constant coefficients. They provide 
also the geometrical extension of Hamiltonian formalism of classical mechanics 
to wave mechanics (see [10]). In applications to physics, it is often called the ray 
theory (see [16]). The same tools apply to the study of the semi-classical limit of 
quantum mechanics and to the high frequency limit of wave equations (acoustic, 
electromagnetic or seismic waves). 

There is a small parameter e > in the theory which is the Planck "con- 
stant" h in quantum mechanics and the wave length or more precisely the dimen- 
sionless ratio between the wave length and the size of the propagation domain for 
wave equations. Most results are only valid in the limit e — )■ 0, but, for simplicity, 
the reader can think of e as a fixed, small enough, number. 

1.1 TOO's 

A pseudo-differential operator (\1/D0) on Mf^ is a linear operator on functions 
/ : M'' — )■ C, := Opg(a), defined using a suitable function defined on the phase 
space, a{x, : © M| ^ C, called the symbol of A^, by the formula (Weyl 
quantization) 

Mf){x) = 7^ / e^^^-y^-^a (^,e^) f{y)dyd^ . 

The function a is assumed to be smooth and homogeneous near infinity in ^. The 
Schwartz kernel^] [/l£](x, y) of is located near the diagonal x = y and is of the 

"'^The "Schwartz kernel" of a linear operator A is the "continuous matrix" of A, we will 
denote it by [A]{x,y) and it is characterized by Af{x) = J^[A]{x,y)f{y)dy. 
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form 

where k{x, z) is a smooth function outside 2 = going to as 2 — oo. 
Simple examples are 

• Op^(l) = Id by the Fourier inversion formula 

• OvAQ = li- 

• Opg(xj) is the multiplication by Xj 

• If X is a positive function with bounded support, the operator Op^{x{0) is 
a frequency filter 

• OPfidCP + = — + V{x): the Schrodinger operator 

• Opg(n(x)|^p) = — e^div grad): the acoustic wave operator. 
The main properties are the following ones which hold as e — 0: 

• Composition: 



Brackets: 



where 



Op, (a) o Op, (6) ^ Op,(a6) 
[Op,(a),Op,(6)]^-Op,{a,6} 
, , ^ / da dh da dh 



is the Poisson bracket. This last property is very important because it 
relates the algebra of \E'DO's to the geometry of the phase space given by 
the Poisson bracket. 

1.2 Wigner functions 

Wigner functions define the localization of energy in the phase space x 
for a wave function u = u{x). They involve the scale e. The Wigner function 
W^{x,^) of u is the function on the phase space defined by the identities 

Va G (M^'^) , / a{x, OW^i^, Odxd^ = {Op^{a)u\u) , 

where {u\v) = J u{x)v{x)dx, or 



I have 

where J^^u{^) = is the e— Fourier transform of u given by 

This means that the marginals of the Wigner measure W^{x, ^)dxd^ are \u{x)\'^dx 
and \F,u{i)\^di. 

1.3 Hamiltonian dynamics and ray method 

Let us consider the wave equation Uu — Cu = where C is an elhptic \1/D0 hke 
the acoustic operator C = div(?i grad). The symbol, usually called the dispersion 
relation, of this equation is u"^ — n(x)||.^|p = 0. To this relation is associated a 
dynamics called the ray dynamics given by the Hamilton equations: 

dxj dH d^j dH 
dt d^j ' dt dxj 

with H = v^ll'Cll- The main result (Theorem 14.11 below) uses the "Hamiltonian 
flow" $t(a;,.^) is the value at time t of the previous differential system ([1]) 
with data (x, ^) at the time t = 0. In the case of an homogeneous medium, 
n = riQ = const ant, I have 

$i(x,o = (x+tv^e/iieiu)- 

In the ray theory, this correspond to the group velocity of waves ^/no. 

The mathematical theory of rays is called the theory of Fourier Integral Op- 
erators and has been developed in the seventies by Hormander and Duistermaat 
(see [TU]) following some pioneering work of Lax and Maslov. A presentation more 
adapted to physicists is given in [16]. Unfortunately, the geometric background 
is rather sophisticated and cannot be presented in a few pages. However, explicit 
formulas in terms of oscillatory functions and oscillatory integrals are available. 

In what follows, I will use the fact that the Green's function G{t, x, y) of wave 
equation admits, in the semi-classical regime (short wave-length), a decomposi- 
tion as a sum of contributions of rays 7 going from y to x in time t: G = ^-y- 

2 Random fields: power spectra and correla- 
tions 

Let / = f{x), X G M'^, be a random complex- valued field with zero mean value. 
Let us denote by E the expectation or ensemble average. 
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Definition 2.1 The correlation of the random field f is the 2-points function 
given by 



The power spectrum of the random field f is the function on the phase 
space given by the expectation of the Wigner functions 



The power spectrum and the correlation contain the same information: 

• The correlation C{x,y) is {{27ieY times) the operator kernel of Op^{p) or 



• is ((27r£:) times) the symbol of the operator whose integral kernel is C. 
Example 2.1 ; the white noise 



Example 2.2 ; a stationary noise on M with e = 1, C{s,t) = F{s — t) and 
Pi{s,u) is the Fourier transform J^{F){oj). 

3 A mathematical model 

I will now discuss the basic mathematical model: it consists of 2 parts : 

• A deterministic wave equation which could be the elastic wave equation or 
more simply here the acoustic wave equation. Because the source of noise 
will be permanent, some attenuation in the equation is needed. 

• A source noise assumed to be stationary and ergodic in time. In seismology, 
this source is usually created by the interaction of the fluids surrounding 
the earth crust (atmosphere or ocean) with the crust itself. This source 
is modeled by a random field which I put on the right-handside of the 
equation. 

For simplicity, I will discuss only the case of a scalar acoustic wave equation 
on some domain in with a random source field / = f{x,t) {t is the time): 



C{x,y):=E{f{x)f{y)) 



Pe := nW'f) . 




C = 5{x-y),Pe = l/{2TieY. 



utt + a{x)ut - Cu = f 



(2) 



where 
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• The field u = u{x, t) is scalar 

• a, the attenuation, is a smooth > function. I will assume for simplicity 
that a is time independent, but it is not really necessary 

• £ is a self-adjoint pseudo-differential operator of symbol —e~'^ll{x,^). Usu- 
ally, Iq is homogeneous of degree 1 which makes C independent of e. This 
will not be the case for dispersive waves like surface waves. Typical exam- 
ples are the Laplace-Beltrami operator of a Riemannian metric on X with 
^o(^)0 = \/ 9^K^)^i^j the acoustic wave operator div {n{x) grad)) with 
hi^ii) = v^n(x)||^||. I introduce Lq := OY<ei}o) = S\/—C 

• f = f{x,t) is a stationary and ergodic (in time) random field with cor- 
relation E(/(s, x)/(s', ?/)) = 6{s — s')T{x,y) and power spectrum p{x,^); I 
assume that p{x,^) has bounded support and that / is real valued and hence 
that p{x,^) is even w.r. to ^. I assume that p is independent of e, this implies 
that the correlation is e— dependent: in particular, T{x,y) « \x — y\/e. 
The source noise decorrelates rapidly as \x — y\ » e. 

The Green's function is the integral kernel G giving the causal solution of 
Equation ([2]) in terms of /: 



Lemma 3.1 The following relation holds: Ga,b{—t) = Cb,a{t)- 

Hence I can (and will!) restrict ourselves to r > 0. 

Using the fact that the source noise is ergodic and stationary, I get the fol- 
lowing result 

Theorem 3.1 The field correlation is given by Equation ^ only in terms of the 
Green's function G and the correlation T of the source noise 



All the work is now concentrated to get a more explicit and more geometric 
expression: this will be done using an expression of the Green's function as a 
sum over rays going from i? to A in time r and using the power spectrum p of f 
which is a semi-classical expression of the correlation of the source noise. 




Our goal is to compute the correlation 






(3) 
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4 The main formula 

Let us denote by Q± (t) the "one-parameter groups" of linear operators generated 
by ±iLo — ea/2: fl+(t)uo is the solution of the differential equation it = (-Lq — 
a/2)u with m(0) = Uq, and similarly for The use of ^l±{t) is a way to split 

the Green function of the wave equation usually given by some "sinus" function 
into 2 exponentials: this way, I reduce the wave equation from an equation with 
of second order in time to a diagonal system of first order in time. 

I will express the result in terms of operators instead of expressing them in terms 
of their kernels (matrices). This gives a much more compact expression! The 
symbol o means the composition of operators while C(r) is the operator whose 
integral kernel (matrix) is Ca,b(t): 



(C(r)n)(A) = / CA,Bi^)u{B)dB . 
Jx 

The main result is 
Theorem 4.1 The correlation is given, for r > 0, as e goes to 0, by 

C'(r) = [fi+(r) + fi_(r)]on, (4) 

with n = Op^{7f) and 

and if a = ao is constant 

I will compare our result (Equations (HI) and ([5])) to the Green's function. 
In the semi-classical regime, i.e. as e — )■ 0, I have 

G{t,A,B) - ^ [(fi+(t) -fi„(t)) oLo^] {A,B) , 

Let us now compute the r— derivative of Ca,b{t): 

-^C(r) - -- (fi+(r) - fi_(r)) o Lq o H . 

In the case of white noise and constant attenuation ao, I know (see for example 
[S] Section 5.1 for a derivation)) that 

which is consistent with the previous semi-classical formula. 
I can now give a more concrete formula: 
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Corollary 4.1 Writing G{t, A, B) as a sum of contributions of rays 7(5) 

with^iQ) = {B,^b) and^ir) = $.(5,^^) = (A^U), I get 

-^CiT,A,B) = J2M,G, , 

7 

with 

1 rO ^ 

2 j-00 

In the case of the white noise p = 1 and a = oq, I recover the formula 

= -l/2ao . (7) 

Let us also remark that, if there is an unique trajectory from 5 to A in time 
r, the prefactor applies to the Green's function itself. It is the case, if I work 
with wave equations with constant coefficients in M". 

The previous formula is consistent with the observations of the paper [22]: 
the correlation Ca,s(t) is not always an even function of r as it is if the source 
is a white noise. The evenness is valid only up to scaling of Ga,b{''')'- 

CsAr) = CaA-^) kGA,B{r) ■ 

The factor k is the ratio of the integrals giving for the ray 7(t) going from B 
to A and 7(— t) going from A to B. 



5 Time scales 

As I see from the general expression of the correlation given in Equation ([3]), the 
proof of the main theorem 14.11 involves the knowledge of the Green's function 
at large times. This is a well known difficulty and the semi-classical expansions 
of the Green's functions are valid up to to the so-called Lyapounov time which 
involves the Lyapounov exponent measuring the rate of instability of the ray 
dynamics. Roughly speaking, the Lyapounov exponent is the smallest number A 
so that the distance between any to rays 71 (t) and 72 (^) satisfies the estimates 

rf(7i(t),72(t))<Ce^*t^(7i(0),72(0) 

with C independent of 71 (0) and 72(0). There is an associated time scale TLyap = 
1/A. On the other hand there is an attenuation time scale for the wave dynamics 
expressed in terms of the decay of the Green's function 

\G{t,x,y)\ <Ce-^/^-' . 

Tatt satisfies the estimate Tatt > 2/infa. The approximation given in Theorem 
14.11 is better when Tatt >> ^Lyap- In particular, this condition is necessary in 
order to get point- wise convergence (i.e. convergence for A and B fixed). 
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6 The use of surface waves for passive imaging 



A remarkable application of the previous tool is to the imaging of the earth crust 
p!5l [25| l26| El Eni El]. This is done using the part of the Green's function associ- 
ated to the surface waves: the earth crust acts as a wave guide on elastic waves 
and these waves follow an effective wave equation. The effective Hamiltonian is 
described now: let us start with the acoustic wave equation uu — diy{n gradw) = 
with the function n coming from a stratified medium n = n(x, z) (here z = 
is the surface) where n is weakly dependent of x (this can be formalized as 
n(x, z) = N{e:x., z) with smooth and e small). Using the adiabatic separation 
of variables u ~ U{ex, z)e*^^'^^ with U weakly dependent of x, I can operate as if 
n was independent of x and I get the reduced equation 

f/« + Opi(A(x,O)f/ = 

where A(x, ,^) is an eigenvalue of the Sturm-Liouville operator 

with appropriate boundary conditions at z = 0. 

From the correlation, I get the ray dynamics of the surface waves and hence the 
effective Hamiltonians A(x, ^). The inverse problem to be solved is the following 
inverse spectral problem: from the fundamental mode (or any other available 
mode) of i^xd in some range of wave numbers |^|, recover n{xo,z). This is the 
kind of well posed inverse problem for which analytical/numerical method can 
be used (see |1]). 



7 A formula for the scattering of random plane 
waves 

I have seen an exact formula for the correlation of the wave field when the at- 
tenuation a is constant and the source noise is a white noise. I will see another 
exact formula in the context of wave scattering by a perturbation sitting in a 
bounded domain of (see [6]). This formula is very general and applies in all 
situations of wave scattering (scalar or elastic waves), i.e. for any medium which 
is homogeneous near infinity: non-homogeneity's lies at finite distances or there 
is a scattering by a bounded obstacle. This calculus was motivated by the result 
of [19], showing that this result is completely general. 

Let us consider for example an acoustic wave equation ([2]) with n = uq outside 
a bounded set of W^. I will consider scattering solutions of the stationary wave 
equation 

div(n gradn) — u'^u = (8) 
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which are of the following form: let us define, for k G M'^, the plane wave 

eo(x,k)=e*^-^. 

I am looking for solutions 

e{x, k) = eo(x, k) + e'^(x, k) 

of equation in M*^, with nok"^ = cj^(d), where e*, the scattered wave, satisfies 
the so-called Sommerfeld radiation condition: 

^•(-.k) = pS7i(e~(^,k)+0(^))..^oo. 

The complex function e°°{x, k) is usually called the scattering amplitude and is a 
signature of the inhomogeneities. The functions e{x, k) are deformed exponentials 
and allow to write an explicit spectral decomposition of our wave operator, which 
is a "deformation" of the Fourier transform. 

Let us look at e(x, k) as a random wave with k = u/^/riQ fixed. The point- 
point correlation of such a random wave C^'^^^{x, y) is given by: 

Cr"(^, y)= f e{x, k)'^^daCk). 



It is proved in [6], section 8, that 

z 

where G{u,x,y) is the stationary Green's function, i.e. the Schwartz kernel of 
(o;^ -|- div(n grad)) . 



8 Conclusions 

I hope to have convinced the reader, even if he is not very much involved in 
mathematics, that it is possible to derive rather explicit asymptotic formulas for 
the correlation Ca^b{t) of seismic noise. The main conclusion is that, in the semi- 
classical regime, even if the source noise is not homogeneous, the field correlation 
is very close to the Green's function; in many cases, there is only a prefactor 
which I computed and which introduces no phase shift. This prefactor vanishes 
if the support of the source noise does not meet the rays from B to A. 
Many other ideas and applications remains to be exploited: 
Is it possible to use the previous tools in order to get informations on the 
source noise? Can I extend the previous calculus to the case where the source 
noise is located on a surface? Can I do something similarly in other regimes of 
propagation, in particular in non-smooth media? Can I get applications of the 
general formula to monitoring? 

^As often, I denote k := |k| and k := k/fc 
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